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Abstract. We study derivatively coupled fermions in axion-driven inflation, specifically m^cj) 2 
and monodromy inflation, and calculate particle production during the inflationary epoch and 
the post-inflationary axion oscillations. During inflation, the rolling axion acts as an effective 
chemical potential for helicity which biases the gravitational production of one fermion helicity 
over the other. This mechanism allows for efficient gravitational production of heavy fermion 
states that would otherwise be highly suppressed. Following inflation, the axion oscillates 
and fermions with both helicities are produced as the effective frequency of the fermion held 
changes non-adiabatically. For certain values of the fermion mass and axion-fermion coupling 
strength, the two helicity states are produced asymmetrically, resulting in unequal number- 
densities of left- and right-helicity fermions. 
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1 Introduction 

In this work we study the behavior of fermionic degrees of freedom during and immediately 
following inflation driven by an axionic field [1,2] (for a recent review of axion inflation see [3] 
and references therein). Specifically, we consider Majorana fermions coupled derivatively to 
the inflaton field via a dimension-five operator. When the fermions have degenerate masses 
and are paired to make Dirac fermions, this is the coupling to the axial-vector current, a 
conserved current in the massless limit (in the absence of anomalous gauge couplings). 

In local thermal equilibrium, this coupling between the fermion fields and the rolling 
axion acts as a chemical potential for helicity, allowing the system to reduce its free energy 
by converting some fermions of one helicity into the other. The universe is not in thermal 
equilibrium during inflation, however, in this work we demonstrate that this coupling has 
a similar effect on the gravitational production of fermions. During inflation, the exponen¬ 
tial expansion of spacetime produces particles from the vacuum as their wavelength becomes 
comparable to the Hubble length. In the absence of the axion coupling, all spin, or helicity 
states are populated with a probability that depends only on the fermion mass, m^, and the 
Hubble rate, H. This gravitational pair production is most efficient for small masses, with 
the production rate for larger masses being exponentially suppressed by the ratio m^/H. In 
the limit that the fermions are massless, no particle production occurs because the theory is 
conformally equivalent to a Minkowski space theory (see e.g. [4]). We show that the coupling 
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to the rolling axion leads to the biasing of this gravitational fermion-production, which en¬ 
hances the production rate of one of the helicities while suppressing the other. Furthermore, 
this enhancement can partially counter the exponential suppression of heavy masses allowing 
for efficient production of fermion states that would otherwise be heavily suppressed. On the 
other hand, in the limit that the fermion mass vanishes, the axion-fermion coupling reduces 
to a boundary term which has no effect on the equations of motion and the effect vanishes. 

The significant feature arising from considering a pseudo-scalar axion held as the inhaton, 
that we will discuss here, is the nature of the derivative coupling to fermions. The specific 
form of the axion potential appears to play a role in the details, but does not change the 
mechanism and the qualitative results of the paper. It is shown that even the final results 
only differ by an 0(1) factor for a significant parameter range in two different forms of the 
inhaton potential we considered. 

Derivative couplings of vector gauge-held degrees of freedom to rolling pseudo-scalar 
helds [5-11] lead to strong growth of one of the polarizations of the gauge helds during de 
Sitter expansion, but the effect of such couplings on fermion helds has been little studied. 
We demonstrate that in the fermionic case there is a similar effect, however, Pauli blocking 
prevents the occupation numbers from exceeding unity and consequently the growth of the 
fermion held is constrained. 

Following the inhationary epoch, the axion oscillates about the minima of its potential 
which causes the effective frequency of the fermion helds to vary non-adiabatically. These 
oscillations cause both helicities to be produced in turn as the effective momentum of the left- 
and right-helicity fermions vanishes. However, the axion oscillations damp as the universe 
expands which reduces the efficiency of successive particle-production events. In combination 
with the inhationary production, this damping can lead to asymmetric number densities of 
fermions with left and right helicities. Moreover, since neutrino helicity is equal to lepton 
number in the standard model, this effect may be relevant for the generation of the matter- 
antimatter asymmetry in the universe - we pursue this question in a companion paper [12]. 

It is well established that fermions can be produced via Yukawa couplings to oscillating 
scalar helds after inflation [13-17]. Preheating in these models is, in some instances, known to 
produce fermions with masses much larger than the inhaton mass and may be important for 
leptogenesis scenarios where the baryon asymmetry is produced by the decay of right-handed 
neutrinos with masses near the grand unified energy (GUT) scale [18, 19]. Further, Ref. [20], 
studied neutrino production during the relaxation and oscillations of the Higgs condensate 
following inflation as a mechanism for the generation of a matter-antimatter asymmetry 
via leptogenesis. The cosmological consequences of the gravitational production of heavy 
fermionic states during inflation is well appreciated [21], and gravitational production has been 
considered as a mechanism for producing super-heavy fermions during inflation [22, 23]. The 
behavior of fermions coupled to axions in the post-inhationary universe was first considered 
in the context of spontaneous baryogenesis [24-28] and has more recently been revived for 
leptogenesis scenarios [29, 30]. Our work on fermion production differs from these existing 
works because we consider both the production during the inhationary period due to the 
axion coupling, which has significantly different phenomenology, as well production following 
inflation during the reheating phase. 

We work in natural units where c = h = 1 and denote the reduced Planck mass by 

M" 1 = V8ttG. 
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2 Background and conventions 


Our notation and conventions for the axion, gravitational, and fermion sectors of the theory 
are as follows. We consider an axion inflationary sector minimally coupled to Einstein gravity 


Sinf = 






( 2 . 1 ) 


where V((j>) is a potential capable of supporting an extended period of slow-roll inflation. 
For concreteness we consider the simplest chaotic inflation model with a quadratic potential 
[31, 32] 

V{4>) = , (2.2) 

as well as the potential from the simplest type of axion-monodromy inflation [33, 34] 

V{4>) = M 3 (V<£ 2 + 4>l ~ 0c)- (2.3) 


In addition to the inflaton sector, we consider a set of fermions 1 Xh described by the action 


S { = d A xyf^g 


IXiae^aV 


-J . 'r aap D li Xip ~ 7,mij(XiXja + xUx] a ) + ^-d^xU^ a^XjP 


Ci 

f 


(2.4) 


where rn t j is hermitian, Ci is real, and / is a mass scale associated with the axion. The 
covariant derivative of the spinor fields is 


D/x — d^ T tO^iab 


_a -b 

a , a 


(2.5) 


where the spin connection is 

mb — 2 (dii&bv — € a &b (dj/CcA d\C C y)& ^ . (2.6) 

We work in cosmic time and ignore metric fluctuations as well as inflaton fluctuations. While 
fluctuations of the axion and gravitational fields are likely to be interesting, our primary 
interest in this work is the behavior of the fermion fields in the classical, homogeneous axion 
and gravitational backgrounds. We work in ‘mostly plus’ convention for the metric with 

ds 2 = —dt 2 + a 2 dx 2 = e a ^e b ^ri a b, (2.7) 


where a(t) is the scale factor and e a ^ are vierbeins. 

We rescale our fermion fields by a~ 3 / 2 to absorb the factor of \J—g in the action and turn 
the covariant derivative into a partial derivative. We write Xi = o^^Xi so tlia,t the fermion 
action reads 


Sf = I d 4 .x 


ixU^aV 


M - ImjixfXja + xUx] d ) + jd^xUe\d a «P Xi p 


■ ( 2 . 8 ) 


1 We work with two-component spinors and consider only fields that transform under the left-handed 

representation of the Lorentz group. Right-handed fields can be obtained simply by complex conjugation, 

that is, if Xi is a left-handed field, xl is right-handed. See Appendix A for some key notations used in this 

work. Refer to [35] for an extensive review. 
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For simplicity, we consider only a pair of fermions, i E {1, 2}. Further, we treat separately the 
case where the mass matrix rriij has degenerate eigenvalues, in which case there is an 0(2) 
flavor symmetry and the fields can be paired to make a Dirac fermion. 

Above and in what follows, Greek letters from the middle of the alphabet represent 
four-dimensional space-time indices, while Greek letters from the beginning of the alphabet 
denote spinor indices. Roman letters from the start of the alphabet denote four-dimensional 
Lorentz indices, while Roman indices from the middle of the alphabet denote spatial indices 
as well as flavor indices. Einstein summation convention is used for spacetime, Lorentz and re¬ 
peated flavor indices, unless otherwise noted. Paired ascending dotted indices and descending 
undotted spinor indices will also be summed. 


3 Majorana and Dirac fermions 

There are subtle differences between the cases where the fermions can be combined into a 
single Dirac fermion and can carry charged currents compared with the case where they 
are Majorana fermions. We treat these cases separately in what follows and then consider 
their vacuum states and quasi-particle numbers in expanding Friedmann-Robertson-Walker 
spacetimes. 


3.1 Non-degenerate fermions: Majorana fermions 

First we consider the case where the fermions have non-degenerate masses and there are no 
conserved currents. Without loss of generality, we work in a basis where rriij is diagonal, and 
the fermions have masses rm. Variation of the action at Eqn. (2.8) with respect to y \a yields 
the equations of motion for these Majorana fermions 

iffad^PduXif} + (3.1) 

It is convenient to expand each held into a Fourier basis, 


Xia(x,t) = 


= E 

A=±l ' 


d 3 &: r 


(2n) 


x ia (Kt)a ik e 


A _zk*x 


+ 2h A a(M)aI£e 


t A — zkx 


(3.2) 


where we have introduced creation and annihilation operators, a,^ k and aj k , which satisfy the 
anti-commutation relations 


a jk'> = (2vr) 3 <5 3 (k - k')5 A y<%, (3.3) 

with all other anti-commutators vanishing, as usual. Here and throughout we use A = ±1 to 
denote the spin or helicity states of the fermions. It follows that 


yt. 

A IOL 


( x , t )= 


A=±l ’ 


d 3 fc 

(27t) 3 


x il ( k ’ f)° T ike * RX + 2 hi ( k > t) a ik e 


p -* k x 4- 


A jk-x 


(3.4) 


We quantize the fields by imposing the anti-commutation relations 

{Xia(x,t),7r^(y,t)} = i<5fh 3 (x - y)<%, 


(3.5) 
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where the canonical momenta of the fermion fields are found in the usual way 


dCf 


p _ u ^f_ _ t . ^Oa/3 

'* v,- a jttO , 


(3.6) 


where here and throughout an overdot on a field denotes a derivative with respect to cosmic 
time, i.e. Xip = dtXi/3■ After inserting the mode expansions, the fields are canonically 
quantized provided that the fermion wavefunctions are normalized according to (no sum on 
*) 


[4( k > *K A i( k > + yiaiK k , = <4 


(3.7) 


Substituting the mode expansions, Eqn. (3.2) and (3.4), into the equation of motion, Eqn. 
(3.1), we find the coupled equations for xp and y^ 


i xfp(k,t) + j^a 0a/3 xfp{k,t) =miy^ a (k,t), 

* + 2/4 CM) - =mi®w(M)- 

We work in a basis of helicity eigenspinors, which satisfy 

3 -hi A = A£a, a = ±1, £_a (~k) =^\(k), 

where t A is a phase that satisfies 


-fc 


Writing the spinors as 

4(M) =A A (f)Uk), 4 f “(M) = 4*(f)£ A (k), 

canonical quantization requires that (no sum on i) 


= 1. 


J2[Xi k miC(t) + Y£(t)Y*(t) 

A 

In terms of these wavefunctions, the equations of motion (3.8) and (3.9) become 

ifdt-i^-x] Xl k {t) + =rmYte(t), 


a 

.k 


f 

^ t ( j.\ v"A . 


i[d t + i -A j Yt(t) - j<f>Y£(t ) =mi^(t). 


(3.8) 

(3.9) 

(3.10) 

(3.11) 

(3.12) 

(3.13) 


(3.14) 

(3.15) 


These coupled equations can be decoupled to give two second order ordinary differential 
equations for each wavefunction 


d't + - id t k x (t )) X&(t) =0, ( d'f + u z x {t) + id t k x (t )) Y.£(t) = 0, 


(3.16) 
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where 


u{(t) =k\(t) 2 + m-, 




(3.17) 


We are interested in particle creation, and so we compute the Hamiltonian density 

77 f = 7rJ.Xi/3 - £f = -a _1 *xL^ Q/3 ^Xi/3 + ^yUiXia + xLxl") - y ^xL^Xi/3- (3.18) 

After Fourier transforming, inserting the mode expansions, and some straightforward algebra, 
the Hamiltonian is given by 


« = £ 


d 3 k 

(2ir) 3 


EHk)a^ + Ft(k)Xf 


V A k a 


t A + T) A *(fc)Ar 


A *n A n A 

a k a_ k 


- ^ ac (fc)(27r) 3 5 3 (0) 


(3.19) 


where 


£, A (fc) =MT A *(f)y, A (t) - x£(t)X*(t)) + mi(X&(t)Y£(t) + Y^(t)X^(t)), 

FHk) =kxX£(t)Y*(t) - ^{Y&Wikit) - 

ET(k) =k x Y^{t)Y^t) + ^(X* (t)Y*(t) + Y*(t)x£(t)). (3.20) 

Note that the presence of the phase prevents the off-diagonal terms in the Hamiltonian 
in Eqn. (3.19) from vanishing. The diagonal and off-diagonal parts of the Hamiltonian are 
related by 

(£ a ) 2 + 4|if | 2 = ( ~k 2 + m 2 ) (x A X A * + Y^*) 2 = ’d. (3.21) 

3.2 Degenerate fermions: Dirac fermions 

In the case where the mass matrix has degenerate eigenvalues there is a global internal 0(2) 
flavor symmetry Xi ~Xj • where O is an orthogonal matrix, provided we assume that 
Ci = C, that is, that the axion couples universally to the fermions. There is a conserved 
hermitian Noether current associated with this symmetry 

=i (xL S^X 2 D ~ xY^Xv) ■ (3.22) 

However, the flavor basis we have been using is off-diagonal with respect to this current and 
thus the flavor states are not eigenstates of the charge operator Q = [ d 3 xJ°. It is more 
useful to work in a basis of states of definite charge. We can achieve this by diagonalizing the 
Noether current. Redefining our fields 


<P^Xi + iX2), ri=-^={xi-iX2), 

diagonalizes the current in Eqn. (3.22). In this basis, the action reads 
St = j d 3 xdt iip^a^d^ipp + - m((p a r] a + ip\i^ a ) 


(3.23) 
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and the SO( 2) flavor symmetry (the part continuously connected to the identity) is realized 
in this basis as the 17(1) symmetry (p —> e %e <p and rj —» e _ * 0 r/, where 8 is the angle that defines 
the 50(2) rotation matrix [35]. 

The spinors p> and r/ satisfy the free-field Dirac equations 


ia^d^p + - rrup^ a =0. 


A 

-h(x,t)=E J (2lr )3 


d 3 fc 

(27t) 3 

d 3 k 


rf 61 =0, 


(3.25) 

p^ =0. 


(3.26) 

expand into plane waves 

as 

,jik. ()/,'> ’ k ' 

J 

(3.27) 



(3.28) 


where a\ b\a, and b are creation and annihilation operators that satisfy the anticommutation 
relations 


*k' J — i w k> w k 

with all others vanishing. It follows that 

*x, t) =5:/. d3 * 

A 

4 ( x ^) 


K,4"'}={&k^i"'} = 5AAo 3 (k-k / ), 


(27r) 3 

d 3 k 


(2tt) 


x^k,t)a*e- ikx + y*(k,t)b& kx 
xl\k,t)b^e-^ + yi(k,t)a^ 


(3.29) 

(3.30) 

(3.31) 


Substituting these mode expansions into the equations of motion, we find the same equations 
of motion as in the non-degenerate Majorana case 


i (a 0al3 d 0 + atg(k, t) + y (frcr 0al3 Xp(k,t) =mj/ At “(k, f), (3.32) 

i (yj% f do + ?'^ y At7 (M) - y<M°^y At7 (k, t) =mx$(k,t). (3.33) 

We proceed in the same fashion as the non-degenerate mass case, and work in a basis of 
helicity eigenspinors, as above 

*«(M) =x£(t)Zx( k), y At "(k,t) = Y k * (i)£\(k). (3.34) 


Then, the equations of motion for the wavefunctions of the Dirac fields are identical to Eqn. 
(3.14) derived earlier for each flavor in the non-degenerate system, 


i (dt - i k ~ a) X x (t) + =mY k x *(t), (3.35) 

i (dt + i\ A) Yf(t) - yj>Y x *(t) =mX k (t). (3.36) 
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However, note that for two fermions there is only a single pair of equations, unlike the Majo- 
rana case where there is a pair of equations for each flavor. This is ultimately because Dirac 
particles are related to their antiparticles by charge conjugation; the 17(1) symmetry relates 
the fields to each other. We can construct the Hamiltonian for this theory analogously to 
before 

H, = £/ i0 5 {EA(*.t)(« , X + '> , ^)-£I“(<:.*)(2») 3 '5 3 (O) (3.37) 

A 

+ F x (k, t)Xi X / X _^ X + F * (k _ t)Ai A* a A 6 A k j ; 


where we have defined 


E x (k) =k x (Y x *(t)Y x (t) - X x *(t)X x (t )) + m(X x *(t)Y x *(t) + X x (t)Y x (t)), 

F\{k) =2 ~k x Xt(t)Y x (t) - m(Y x (t)Y x {t) - X x *(t)X x * (t)), 

Er(k) =2 k x Y x *(t)Y x (t) + X x (t)Y x {t) + Y x *(t)X x *(t)), (3.38) 


which is simply the sum over the flavors of the Hamiltonian at Eqn. (3.19) where nii = m. 
The relations between these functions for Majorana and Dirac Hamiltonians are 

E^™ c {k) = 7if ajorana (A;), (3.39) 

i? Dir ac ( /c ) = 2i? Majorana^^ ( 3 . 40 ) 

and the corresponding consistency relation in the Dirac case becomes 

{E x f + |F A | 2 = ( k 2 x + m 2 ) (X&X£ + Y X Y X *) 2 = (3.41) 


3.3 Vacuum states 

The Hamiltonian in both Eqn. (3.19) and Eqn. (3.37) can be diagonalized at any given time to 
by choosing the wavefunctions X k and Y k so that the function F k vanishes. A straightforward 
calculation shows that F k = 0 whenever 


oc 




JL h 


OC 




(3.42) 


where (j) is an arbitrary phase. Here, we are interested in modes that start in the Bunch- 
Davies vacuum state during inflation. In order to determine the initial condition, we look 
for Wentzel-Kramers-Brillouin (WKB) solutions to the equations of motion. Assuming that 
|(h/(u 2 | <C 1, it is straightforward to derive the lowest order WKB solution to Eqs. (3.14) - 
(3.16) 


X, 


Ut) = Aa \ 1 + + B iX J 1 - 




w A 


+ A J 1 + ^ 


k\ _ 


W A 


WA 
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Y X *(t) = - A A 


i f u)\dt 


(3.43) 









where as above 


w\(t) =k\(t) 2 + rnf , k x (t) = + j4> \ . (3.44) 

Substituting these solutions into Eqn. (3.20) we obtain 

Ehk) - |A a | 2 ), F?(k) = 2uA* x B* x , 

Er(k) =h(\A iX \ 2 + \B iX \ 2 ) + u x (\Bix\ 2 -\A iX \ 2 ), (3.45) 

for the Majorana case, while substituting them into Eqn. (3.38) gives 

Ef(k) =2u x (\B lX \ 2 - |M lA | 2 ), F^k) = 4u ,A* X B* X , 

Er (k) =2k x {\A iX \ 2 + \B lX \ 2 ) + 2co x (\B iX \ 2 - |A a | 2 ), (3.46) 


for the Dirac case. On the other hand, canonical quantization requires that the coefficients 
satisfy 

£ [lUiAl 2 + |Ba| 2 ] = \. (3.47) 

A 


Thus, in order for the fields to be canonically normalized and in the lowest energy state (the 
state with the lowest vacuum energy density), we are required to take 


Ai\ — 


1 

2 ’ 


B iX = 0 . 


(3.48) 


Notice that in the absence of the coupling to the axion (when k\ = kX/a), the first term 
in the vacuum energy in Eqn. (3.45) - (3.46) vanishes when summed over spins. When the 
coupling to the axion is turned on, this term no longer vanishes, and leads to a splitting of 
the vacuum levels for the A = +1 and A = —1 helicity states. 

The choice of initial conditions at Eqn. (3.48) diagonalizes the Hamiltonian at some 
initial time, but in general the equations of motion drive F^{k) away from zero at times 
t > to. At these later times, the Hamiltonian can be re-diagonalized after performing a 
(time-dependent) Bogoliubov transformation on the creation and annihilation operators. 

We begin by constructing the Bogoliubov transformation for the Majorana case, by first 
writing the Hamiltonian in the following form 



( E*/2 ( a t \ 

I-BaMW ~E X /2 ) W.J 


This Hamiltonian can be diagonalized by choosing 

«k =«fc(i)ak + /3k{t)t- k a ^ k , 


£X ac (27r) 3 <5 3 (0) 


(3.49) 


(3.50) 


where ak(t) and /3fc(i) are complex functions of time that depend on the magnitude of the 
wavenumber. Under this Bogoliubov transformation, the creation-annihilation operator pairs 
become 


a\a k = \a\ 2 a\a k + |/?| 2 (t-fci* k )a- k a_ k + a*/3L- k a' k a_ k + /3*at*_ k a_ k a k , 

= ( a *) 2a U a l + (P*) 2 (E_ k t* k )a k a_ k + a* fi* t,*_ k ai_ k ai_ k a_ k + /3*a* F k a k a\, 
a k a_ k = a 2 a k a_ k + /3 2 (i k t- k )a^_ k al + a(h k a k a\ + f3cu- k a ] _ k a_ k . (3.51) 
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The off-diagonal part proportional to a k a_ k is 


-E/3*oa*_ k + F(P*)\i*_ k i* k ) + F*a - 2 , 


(3.52) 


where we supressed helicity and flavor indices. Diagonalizing the Hamiltonian amounts to 
setting the off-diagonal terms equal to zero. Solving for a is now straightforward, setting 
Eqn. (3.52) to zero yields 


where we used i_ k 
leads to 


a = 




(3.53) 


— Lk, and took E 2 + 4|F| 2 = tc 2 /4. Using the relation |a| 2 + \/3 \ 2 = 1 


8|F| 2 

u:(2E + u) 




(3.54) 


where we used the fact that |^4u| 2 + I^aI 2 = 1/4 when not summed over spins. 

We can proceed in a similar fashion for the Dirac case, where the Hamiltonian can be 
written as 



Fi c (2vr) 3 <5 3 (0)| , 

(3.55) 


where 


#2x2 


( E x /2 -F x (k) Xi x \ 
\-F*(k)\ l x * -E x /2 )■ 


(3.56) 


Note that the Dirac Hamiltonian is comprised of two identical copies of the Majorana 
Hamiltonian, as expected. The form of the Bogoliubov transformation can be read off as 


Ok =&k(t)a£ + 


b-k = ~ Pk(t) a t + a k (t) b l 


k > 


(3.57) 

(3.58) 


as, for example, in Ref. [15]. The complex phases t k are absorbed into the parameter /3 in this 
case, since - in contrast to the Majorana case - we do not need to invert the wavenumber of 
the operators as defined in the Bogoliubov transformation, so the antisymmetric nature of i k 
for the relabelling k —> — k is irrelevant. 

Performing the Bogoliubov transformations on each part of the Dirac Hamiltonian, we 
obtain (suppressing flavor and helicity indices for clarity) 

a\a k = \a\ 2 a\a k + \/3\ 2 b_ k tf_ k + a*Pa\.b[_ k + a/3*b_ k a k , 

^_ k b- k = \/3\ 2 a k a k + | a\ 2 tf_ k b_ k - P*aa k b~ k - f3a*b^ k a\, 
b- k d k = -aj3a\a k + a(3b_ k tf_ k - /3 2 aj.i>l fc + a 2 b_ k a k , 

a[bl k = -a*P*a\a k + P*a*b- k b ] _ k + (, oc*fa\b ] _ k - {/3*) 2 b. k a k . (3.59) 

The off-diagonal terms (proportional to b_ k a k ) in this expression are 

(2 EaP* + Fa 2 - F*(/3*) 2 )b_ k a k . (3.60) 
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Diagonalizing the Hamiltonian requires that these terms vanish and gives 

a= - (E + „/ 2 y 

Using the relation \a\ 2 + |/3| 2 = 1 we get 

\/3\ 2 = 4|H| 2 , 

as in the Majorana case. 

We use these time-dependent operators (the a) and b') to dehne (time-dependent) Fock 
spaces built from the zero (quasi)particle states 

a^(t)|0,t) = &^(t)|0,t) = 0. (3.63) 

As noted above, a state that starts off at time to as the vacuum state, i.e. a state containing no 
particles, evolves to a state which at a generic later time contains particles. This is because the 
particle number at a given moment is provided by operators which instantaneously diagonalize 
the Hamiltonian. This instantaneous particle number can be computed by projecting the full 
mode-functions onto an instantaneous WKB basis (Eqn. (3.43)), and extracting the negative 
frequency parts. 

For both Majorana and Dirac Fermions, the particle number can be evaluated as 

n k = \Pk\ 2 = MB£\ 2 . (3.64) 


(3.61) 

(3.62) 


We express this quasi-particle number in a more useful form, by noting that the particle num¬ 
ber can be extracted by comparing a given solution of the Dirac equation to an instantaneous 
WKB solution of the form Eqn. (3.43). The square of the coefficient of the negative frequency 
mode is the particle number and can be extracted via 


n ik = 


1 


u\(k\ + u\) L 


I Y^ I 


+ u 2 \X? k \-2 u&iX&Xg) 


(3.65) 


Note, however, that this particle number is only well defined where the WKB approximation 
is valid. That is, during times where |ch/u;?| <C 1. 


4 Inflationary solutions 


During the inflationary epoch, the accelerated expansion of spacetime leads to particle pro¬ 
duction of all fields that are not conformally invariant [4], In this section we derive analytic 
solutions for the quasi-particle number by directly solving the Dirac equation for the fermion 
wavefunctions. 

During slow-roll inflation we can treat the inflationary spacetime as de Sitter space to 
a good approximation. In this section, we work with conformal time, which is taken to be a 
negative increasing quantity during inflation 


dT 



d t 
a 


1 

aH' 


(4.1) 


We assume that the axion rolls at an approximately constant rate in cosmic time during 
inflation (j) “C and we introduce the parameter 


■Q 


7 H 


const. 


(4.2) 
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In conformal time, the equation of motion, Eqn. (3.16), for the spinor fields becomes 


##k +=(i +aw «* + (*»+ * + *i3L + #»i)«£_a, 


ikX 


s± 


— T 


1 If mj 


4 t 2 t 2 [H 2 


s± 


(4.3) 


where we have treated H ~ constant and 0 ~ constant, dropping time derivatives of these 
quantities, and we have rescaled 


v A 

=- * 


a 


(4.4) 


Redefining the time coordinate to u = 2ikr puts the equation into standard form of the 
Whittaker equation 


^uAk ( X + ) Ak. + ( — T + -T H 77 ( Yvk + "i? 


A (1 


u V 2 


1 1 


1 / m 


4 4u 2 u 2 Vif 2 


Ak = o, 


with solution 


where 


Ak = A ik W *A 2ikT ) + B ik W -^A - 2ikr )> 


K = “ A (n+^)> A* = - ( ^ # 




(4.5) 


(4.6) 


(4.7) 


We make use of the fact that in the limit z —>• oo, the Whittaker functions have limiting forms 
[36] 

W^{2iz) ~ e~ iz (2iz) K , (4.8) 

while in the same limit, the WKB solutions to the Dirac equation, Eqs. (3.43), are 


lim 

kr — y —oo 




'1--e 

w+ 

k- 


y/2k 


lim W1 - — e -iS“- dt =V2e ie2 e^{2ikTAe~ ikT , 


kT — y —oo 


UJ- 


(4.9) 


(4.10) 


where and 02 are irrelevant phases. Thus, choosing the Bunch-Davies vacuum implies we 
should set B^ k = 0 in Eqn. (4.6). 

Next, we find YA* in terms of XA and thus A k by using the Dirac equation, Eqn. (3.14). 
In terms of 4>^ k = Y ik */^/a and with u = 2ikr, Eqn. (3.14) is 


A k = - 


iu 


mt/H 




2 u 


mi/H 


d u + A 


k 1 

u 2 


Ak- 


(4.11) 


To solve this equation, we make use of the fact that the Whittaker functions satisfy the 
identities [36] 


d 

Z dz Z 


e ^ z z -K-U 


W k Az))=(1+^k) ea 2 z K W K -n :fl (z), ( 4 - 12 ) 

(e-**z K ~ 1 W K Azj) =(-l) n e-i 2 ^+"- 1 lT K+ri , /i (z), (4.13) 
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where (o) n is Pochhammer’s symbol, which satisfies 


(a)o = 1, (a) n = a(a + l)(a + 2)... (a + n - 1) = 


T(a + to) 

r(o) ' 


(4.14) 


For n = 1, Eqs. (4.12) and (4.13) are 


* (-<9* - \ + W Kili (z) = - f Q - k) - A\ W k -x 

z (~dz + \~ W *.A Z ) =W»s+l,/xW» 


(4.15) 

(4.16) 


and note that, for A = — 1, 


k - n = 


rrq 

i7 2 ' 


(4.17) 


The full solutions to the Dirac equation that match onto the vacuum in the UV are then 
given by 


X iC=r) = ~^W_, 2 _ Ull (2ikT), X~ k (kr)= A “ 


Y ,t'( kT ) = 


^ W i+‘«J 2ikT '>’ < 4A8) 

-1 +H>,A 2ik A- 




It remains to hx the overall normalization of these modes according to (no sum on i) 




= i. 


(4.19) 


Making use of Eqs. (4.12) and (4.13), the fact that the Whittaker functions satisfy the con¬ 
nection formula 


W K ^{z) = W K ,^(z), 

and have Wronskian, 

W[W K ^(z), W- Ki p{e**z)} = e*™ K , 
the normalization condition at Eqn. (4.19) reduces to the condition 

lAtM” = Mil.-** =1- 

Therefore, our canonically normalized fermion wavefunction solutions are 


I , . imie l6 e ^ 

X tk( kT ) = ~WT fnTTZ 


H y/2 kr 


(2 ikr), X.Akr) = 


e* e e2^ 

V^r Wl+i ^ 


jf)’ ——. jQf —19 

p o'-' I'm • c> LXJ p o u 

5 -* (tT) - 5 ‘ r(fcT) --TW W - • 


(4.20) 

(4.21) 

(4.22) 

(2*fcr), (4.23) 

r !+*>,„ (2ifcr), 


- 13 
















m/H (3 


Figure 1 . We display the dependence of the quasi-particle number (for modes with k < aH) on the 
mass of the fermion (in units of the hubble rate) and the coupling strength to the rolling axiom In 
the left panel we show the particle number for fixed values of i? = 10, 0.05,10“ 3 (curves from from 
right to left) as the fermion mass is varied. In the right plot we show the particle number for fixed 
m/H = 0.01,0.1,0.5 (curves from top to bottom) as the coupling to the rolling axion is varied. In 
both cases, we have chosen d > 0 so that n// (red curves) is enhanced while n~j/ (blue curves) is 
suppressed. 


where 0 and 6' are arbitrary phases. 

In the absence of the axion coupling, there are well-known solutions to the Dirac equation 
in de Sitter space, see for example Refs. [21, 37-39]. These results can be obtained after some 
algebra from Eqn. (4.23) by making use of the relation between the Whittaker and Hankel 
functions [36] 


W 0 ,„(2z) = 

and the recursion relations for the Whittaker functions 


(4.24) 


2hW k ^{z) - VzW K+ i ^ + i(z) + VzW k+ ^_i 2 (z) =0, 
( K - M 0 yfzW K _ | )M+ i (z) + 2 nW^z) - + y, - 0 VzW K _i jtl _i(z) = 0 . 


(4.25) 

(4.26) 


Using the analytic expressions in Eqn. (4.23), we derive analytic expressions for the 
quasi-particle number during inflation. We are interested in modes late in inflation that have 
left the horizon. For modes which satisfy k/aH <C 1, and assuming ^ 0, we find 2 


n 


± 

ik 





sinlr 

71 (\/w + ± 

sinh 

2 '(VS+tf 2 ) 


(4.27) 


2 When |i?| rrii/H, Eqn. (4.27) is an excellent approximation to the particle number for modes that 
satisfy k/aH < |#|. 
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Eqs. (4.23) and (4.27) are the main result of this section. 

In the absence of the coupling to the axion, production of both fermion helicities is 
symmetric, as expected, and highly suppressed for fermions with masses larger than the 
Hubble rate. For small masses, the occupation number approaches its maximum value of 
1/2 as rrii/H —> 0. However, note that for m, = 0, the theory is conformally equivalent 
to Minkowski space, and no particle production occurs. 3 When the coupling to the axion is 
switched on, the particle production here is asymmetric between the helicity states. For d > 0 
(t9 < 0), particle production of the A = + (A = —) helicity state is enhanced while particle 
production of the A = — (A = +) mode is suppressed. In Fig. 1 we display the dependence 
of the occupation probability on the quantity t?, defined in Eqn. (4.2), which encodes the 
effective axion-fermion coupling, and the mass of the fermion in units of the Hubble rate. 

In the case where the fermions are degenerate in mass and combined to make a single 
Dirac fermion, we note that while the axion coupling leads to the asymmetric production 
of helicity states, this coupling keeps the number of particles equal to the number of anti¬ 
particles, so that the overall Noether charge is conserved. In this case, if left-handed particles 
are produced, so are equal numbers of right-handed anti-particles. 

5 Post-inflationary evolution and reheating in quadratic chaotic inflation 

So far, we have only solved for long-wavelength modes of the fermion field in a basis of quasi¬ 
particles states. These states do not coincide with an observable basis until the Universe is in a 
matter- or radiation-dominated phase. 4 In the absence of instantaneous reheating, at the end 
of inflation the inflaton oscillates about the minima of its potential. As the axion oscillates, 
the effective frequency of the fermion modes (see Eqn. (3.17)) changes non-adiabatically when 
the effective wavenumber, k\(t), vanishes 

k\{t) = -A + = 0. (5.1) 

a f 

Around these times we expect particle production to occur. Note that for one of the helicity 
states, this non-adiabatic evolution occurs for the first time during inflation and corresponds 
to the particle production considered earlier in Section 4. While exact analytic calculations are 
difficult in this regime, we can gain some insight by making use of the WKB approximation. 
This is the calculation we turn to next. The equations of motion for the Majorana and Dirac 
cases are identical, therefore we treat them uniformly, unless otherwise noted. 

In order to simplify notation for the rest of this section, the comoving wavenumber k and 
the fermion mass m* are measured in units of the axion mass m</>. Similarly, time is measured 
in units of m/) 1 and / and <j) are measured in units of Mp\. We also drop the subscript i on 
the fermion mass, because we focus on only one flavor of fermion. 

5.1 Static universe approximation 

In order to gain intuition about the post-inflationary processes in this model, we first work 
in the limit where the universe does not expand. In this static universe limit, the oscillations 

3 Eqn. (4.27) is derived using in the inflationary solutions to the Dirac equation (Eqn. (4.23)) in Eqn. (3.65) 
and taking the limit fcr —» 0. This limit does not commute with the limit m; —> 0. Taking rm —¥ 0 and then 
taking the limit kr —> 0, one finds the result n± k = 0, as expected. 

4 In the case where the axion coupling is absent, the particle number produced during inflation coincides 
with the particle number evaluated later during the matter-dominated phase that results from the axion 
oscillations about the minima of its potential. 
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of the axion about the minima of its potential are harmonic, cj)(t) = <f>o sin(i), and we can use 
the methods outlined by Peloso and Sorbo [15] to study the evolution of the fermion particle 
number. Peloso and Sorbo used a Yukawa coupling of the inflaton to the fermion sector, 
however, the equations for the axial coupling considered here can be mapped directly onto 
the Yukawa case by identifying 

m(t) I Yukawa k\(t) \ axion ? k | Yukawa t 771 1 axion • (5.2) 


Proceeding analogously to Ref. [15] we look at the points of non-adiabatic behavior, where 
k(t*) = 0. These points occur whenever kX = —(C/f)(f) for the wavenumber of interest. For 
the rest of the static universe analysis we choose A = 1 without loss of generality. We thus 
consider modes for which k < {C/f)<f> o, where = 0 is possible, and work in the limit 

j (/>0 1 in order for the WKB approximation to be valid. Defining the time t* by k(t*) = 0, 

we expand the effective wavenumber near these points as 


k « L, (t - U) = -A(t - £*), (5.3) 

dr 

where A = j(f>o. These zero-crossings happen multiple times and differ based on helicity 
and wavenumber, so we should be writing t* = tf(k). However, we keep t*, to not overly 
complicate our notation. Introducing the parameters 


P = 



m 


(A 2 - fc 2 ) 1 / 4 ’ 


2 = 



(t - U) 


{A 2 - k 2 ) l/A {t - f*), 


the equation of motion for the fermion wavefunctions is approximated near t* as 


+ {p 2 + i + z 2 )Xi k — 0 , 
^Xik + (p 2 ~ i + z 2 )^ = 0 • 


(5.4) 


(5.5) 

(5.6) 


If the fermion held started in the vacuum state initially (t < t*), then the particle spectrum 
after the point where k = 0 for the first time (t > t*) is given by [15] 




m 


dk/dt 


(5.7) 


where the superscript (1) denotes the first k zero-crossing, henceforth also called a production 
event. We call the combination 7r p 2 the fermion production exponent at the i’th production 
event. This is independent of i in the static universe case, but not once we take into account 
the effects of the expansion. In the static universe approximation, the particle number after 
the first production event becomes 

„<" = Ht-'A) ’* <j4 = 7*>’ (5 .8) 

\ 0 ,k>A = j4>o. 

Fermion production is exponentially sensitive to the square of the fermion mass and 
the inverse of the axion-fermion coupling constant. Production is therefore significant when 


- 16 















urn 2 /A < 1 and effectively shuts off when nm 2 /A > 3. A important characteristic of the 
present model is that each wavenumber undergoes particle production at a different time, 
defined by k = 0, which for the static universe case is simply i* = arccos ( — k/A). In the 
Yukawa-coupled model studied in Ref. [15], particle production happens uniformily for all 
wavenumbers k when the effective fermion mass crosses zero. 

We have thus far only considered the first instance of fermion production where many 
modes are uniformly populated up to a maximum wavenumber set by the axion coupling. 
There are many oscillations of the background axion field, and we are primarily interested 
in the occupation numbers after a very long time. For bosons, occupation numbers grow 
exponentially due to Bose-enhancement and parametric resonance. However, fermions obey 
Fermi-Dirac statistics which limits their occupation number to be less than unity. Parametric 
resonance 5 and Pauli blocking in the fermion case leads to chaotic filling and emptying of 
fermion states. In the static universe approach, the time averaged particle number count 
can be evaluated semi-analytically. Following the analysis of Green and Kofman [13, 14], we 
define a locally averaged particle number count 

1 

n k {t) = - J n k (t)dt = F k sm 2 (v k t ), (5.9) 


where T is the period of cj)(t ). The amplitude F k and frequency v k are 

F - - - < 5 - 10 > 


The function is the first fundamental solution of the equation of motion, meaning that 


r(i) 


t =0 


= 1 and dtZ 


(i) 




t =o 


= 0 . 


The WKB method used in Ref. [15] provides a different construction for the particle 
number after one full inflaton oscillation. The two different analytical approximations to the 
fermion number produce a final result of the same form n k ~ F k sm 2 {y k t). However, these 
analytic approximations are originally constructed for different regimes. Peloso and Sorbo 
worked in the broad resonance regime 6 and used the WKB approximation, which produces 
results only for k < {C/ f)<j> o, which is not a restriction for the analysis of Green and Kofman. 
However, the method of Peloso and Sorbo extends to the expanding universe case, as we 
see in later sections. The agreement between the two formulas is excellent for (C/f)<po > 1 
(the “broad resonance” regime) and there is an increasing difference as the coupling decreases 
(the “narrow resonance” regime). We tested these approximations against results obtained by 
solving the Dirac equation numerically on the homogeneous axion background for the broad 
resonance regime, where both formulae give identical results. Fig. 2 shows that the agreement 
between the approximate and numerical results is exact at the points t = 2i rn (where n is an 
integer). That is, the approximation exactly reproduces the particle number after an integer 
number of inflaton oscillations. Fig. 3 shows the evolution of the fermion particle number for 
a system that is set at the threshold between broad and narrow resonance (C/f)4 >o = 1. 


5 It is true that parametric resonance is usually associated with exponential growth which can occur for 
Bose fields. We use the term to describe the fermionic case, acknowledging the significant differences due to 
Pauli blocking, in the spirit of the discussion found, for example, in [14]. 

6 We define the terms broad and narrow resonance in this context based on the magnitude of (C/f)(/)o, 

inspired by the corresponding distinction for bosonic parametric-resonance. 
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Figure 2. Exact (red) and approximate (red) particle number for m = 1, j4> o = 10 and k = 1,10,12 
(top left, top right and bottom left respectively). The blue dots correspond to t being an integer 
multiple of 27r. The values of (black solid) and (blue dotted) are shown in the bottom right 
figure. The three pairs of dots show the values of F\ (red dot) and (green dot) that correspond to 
the three panels of this Figure. 



t 



k 


Figure 3. Left panel: Evolution of the particle number for m = 1, j</>o = 1 and k = 1 (blue), k = 0.7 
(red). The black lines show the particle number envelope given by Eqn. (5.9), (5.10). The vertical 
lines correspond to points where t is an integer multiple of 27r. Right panel: Envelope function F ^ 
(black) for m = 1, j</>o = 1 along with the particle spectra for t = 2tt (blue), t = Fit (red) and t = 207r 
(green). 
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Figure 4. Particle number for A = 50 and m = 1,3,5,7,10 (blue, red, green, black and brown 
respectively). In the nk( 7r) plot the green curve is multiplied by 3, the black curve is multiplied by 
10 and the brown curve is multiplied by 300. In the rik{7 tt) and nk(8ir) plots the brown curves are 
multiplied by 10. All curves are shifted vertically by 1, 2, 3, and 4 for visual clarity. The vertical line 
at k = A = 50 shows the maximum produced wavenumber according to the WKB method. 


In Fig. 4 we show the results of numerically solving for the evolution of the particle 
spectra for different fermion masses with A = C/fcj) o = 50, chosen to be well within the 
region of validity of the WKB results. From comparison of the blue and red curves in the 
upper right panel of Fig. 4 with the corresponding curves in the lower panels, it appears that 
the mass dependence of the final particle spectrum is weaker than that of the initial particle 
spectrum (the particle spectrum after the first zero crossing of the effective wavenumber, 
k = 0). Further, for a range of masses, in our example m < 5, the particle spectra after an 
even number of zero-crossings of k are more similar than their initial spectra would suggest. 
Lighter fermions are more easily produced in the first oscillations, yet the final occupation 
numbers appear to be similar, with the occupation number of several wavelengths reaching 
unity. This effect is due to the nature of fermion production. The occupation number of 
each mode ilk cannot exceed unity, therefore systems with initial particle spectra of nk ~ 1 
can be either reduced or stay almost constant. This behavior causes the particle spectrum 
to develop fine band structure as parametric resonance develops, which can be averaged out 
when calculating the total particle number. 

For light fermions with m < 1, nk ~ 1 for all excited modes after the first production 
event. In these cases, the second zero-crossing can only de-populate modes, since modes that 
are already at nk ~ 1 cannot be further populated. This behavior results in a decrease in 
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Figure 5. Left panel: The average occupation number as a function of the number of zero-crossings 
for A = 50 and m = 1,3, 5, 7 (blue-dotted, green, red and black respectively) in the static universe 
approximation. Right panel: The average occupation number as a function of the fermion mass. 
The blue curve corresponds to the first zero-crossing, while the red and green curves correspond to 
(rik{t = 177t)) and (■ rik(t = 187t)) respectively. 


the averaged particle occupation number after the second production event; modes that were 
rtfc ~ 1 after the first event end up with occupation numbers nf. <C 1 after the second event. 
This effect can be seen in Fig. 4, especially for the blue dotted curve. 

Systems with heavier fermions have an initial particle production rate that is much less 
efficient, resulting in smaller initial particle spectra, rik <C 1. However, for these modes, 
each particle production event is more likely to create more particles than destroy them, 
resulting in an occupation number that grows with time. Given enough time these models 
can produce heavy fermions more effectively (at least for certain wavenumbers) than their 
initial production rate would suggest. In these cases, the final result is a particle spectrum 
that has broadly similar features to its lighter counterparts, as can be seen from the green 
curve in Fig. 4. The important difference between the final particle spectra for fermions of 
different masses is in the range of wavenumbers that reach occupation numbers of ~ 1 
rather then in the amplitude of the occupation number (which is < 1 by definition for 
fermions). 

This picture changes dramatically once the expansion of the universe is taken into ac¬ 
count. We show in the next section that the expansion introduces new effects which provide 
the means for generating a significant left-right helicity asymmetry in heavy fermions. Before 
proceeding to the expanding universe case, we calculate the average occupation number per 
unit occupied volume in momentum space 


(n(t)) 


f 0 °° n k (t)k 2 dk 
/ 0 fcmax k 2 dk 


U3 

A/m! 


nk(t)k 2 dk, 


(5.11) 


where k max is the maximum wavenumber that is excited according to the WKB approxima¬ 
tion. For large values of the coupling, to a good approximation no particle production occurs 
for k > fcmax, and Eqn. (5.11) provides an estimate of the (weighted) average occupation num¬ 
ber of each bin in fc-space. We plot this quantity in Fig. 5 for A = 50. However, we caution 
that this quantity is dependent on fc max itself. For small masses, where particle production is 
more effective, we see a large variation in the particle number after integer and half-integer 
numbers of axion oscillations. However, this is not true for larger masses where, as discussed 
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above, after each axion-crossing the particle number is predominantly rising. In the expand¬ 
ing universe case, this oscillating behavior of the average particle number is suppressed, since 
particle production after each production event is less efficient than the previous one due to 
the fact that the amplitude, and thus velocity of the axion oscillations, is damped by Hubble 
friction. 

We end this section on the static universe approximation with a comment about the 
dependence of the maximum produced wavenumber on the coupling strength. According to 
the WKB approximation, particle production occurs only for k/A < 1. However, from Figs. 
2, 3, and 4 note that this limit is increasingly violated as the coupling gets smaller. For 
couplings approaching unity there is significant particle production for k/A > 1, which is 
absent for larger couplings 1. Parametric resonance allows for the production of these 
high-/;; modes. 


5.2 Expanding universe 

The static universe approximation is a useful simplification that allows us to build intuition 
about the relevant physical processes without the complications of expanding space. How¬ 
ever, away from the individual particle production events, it is not a good description of the 
processes occurring during and after inflation. Following inflation, the oscillation frequency 
of the background inflaton held is typically not much faster than the expansion rate, and con¬ 
sequently it is important to take expansion into account. As we demonstrate in this section, 
there are two main differences between the static- and expanding-universe cases. Both the 
range of wavenumbers that get excited and the particle production strength are progressively 
reduced by the expansion, in contrast to the static universe, where they remain constant. 

The main goal of this section is to estimate the final occupation number of the fermion 
states. Additionally, we calculate the initial and final fermion production asymmetry between 
the two helicity states. Since one helicity is being produced for the first time during inflation, 
while the other is produced only during preheating, we will demonstrate that unequal numbers 
of each helicity are produced in the overall evolution. We show that, because the behavior of 
both helicity states during the preheating phase is comparable and highly stochastic, the early 
behavior determines the ultimate fate of the asymmetry. Provided this initial asymmetry 
is large enough, we demonstrate that, at least in the model at hand, it can be preserved 
throughout the preheating period. This is quantified in Section 5.2.3. 

Analogously to the WKB-based static universe approach, particle production in expand¬ 
ing spacetime occurs in distinct incidents at the times when fc(t*) = 0, corresponding to 

-A-A = -^(f*). (5.12) 

a{ U) / 


Fermion production occurs in a very narrow region around k = 0, therefore the expansion 
of the universe can be neglected during the production event and the static universe result 
of Eqn. (5.7) can be immediately applied in this case. Different wavenumbers are excited at 
different times, according to Eqn. (5.12) and as shown schematically in Fig. 6. The dependence 
of the particle production exponent on helicity and wavelength is hidden in the term 


d k 

C 

d t 

u~ 7 


0(f*)i7(f*) + <£(£*) , 


(5.13) 


where the time t* depends on k through Eqn. (5.12). 
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Despite describing particle production starting from a vacuum state, the form of 
given in Eqn. (5.7) is useful for calculating statistical properties of the particle spectrum, 
even after many production events and the emergence of a fine band structure, as shown in 
Fig. 4 for the static universe case. We generalize this form in the following sections to include 
the case of a(t) ^ const and multiple productions. We consider rik as given in Eqn. (5.7) to 
be an indicator of the strength of any single particle production event. 


5.2.1 During inflation 


As discussed earlier and in Section 4, the axion-fermion coupling is active during the infla¬ 
tionary epoch. For one of the helicities (in our case, with (ft < 0 or equivalently id > 0, the 
A = +1 helicity) k\ may vanish during inflation resulting in particle production during the 
inflationary epoch itself. In this case, the time of the production event, where fc(t*) = 0, is 
found to be 






(5.14) 


where we used the approximation of exact de-Sitter space (H = const.) together with the 
slow-roll approximation eft = 0. We see that for (ft < 0 only the positive helicity state (A = 1) 
can be excited during inflation. 

The resulting fermion number with momentum k can then be evaluated from Eqn. (5.7) 
yielding 


( irm 2 A 

nfe=exp r^J- 


(5.15) 


We compare this result to the exact solution derived for particle production during inflation 
given in Eqn. (4.27). Working in the limit id 3> 1 (actually we do not need it to be much 
bigger than one, simply a few times larger), as well as m = O(H), Eqn. (4.27) becomes Eqn. 

(5.15) for the growing (A = +1) state. This WKB analysis agrees with the exact solution 
obtained for the modes during inflation, in the appropriate limit of large coupling (or id 3> 1), 
where the WKB method is applicable, leading to a unified treatment of the inflationary and 
post-inflationary fermion production for the parameter space of interest. It is worth noting 
that the WKB method employed here is valid only in the large-coupling regime, and therefore 
it cannot capture the particle production occurring solely due to the expanding space-time. 

The mode that is excited at the end of inflation can be read off from Eqn. (5.12) to 
be k = — (C/f)X(ft(t en d)a(t en d)- In the case of chaotic inflation with a quadratic potential 
< ft(j'end)Q j (tend) ~ 0.7. 


5.2.2 After inflation 

Following the inflationary epoch, the axion oscillates about the minima of its potential. The 
particle production due to each zero-crossing of k (neglecting the emerging band structure 
due to stochasticity) is governed by Eqn. (5.7) up to a maximum wavenumber given by 

^max = Ty(a^>)| max , (5.16) 

We consider first the A = — 1 helicity state which, for chaotic inflation with our conven¬ 
tions, is excited only after inflation has ended and the axion has crossed zero. The maximum 
particle number is produced for low wavenumbers, since \dk/dt\t t is a decreasing function of 
k. For low wavenumbers, k <C C/f, the condition k = 0 becomes 0(f*) ~ 0 and particle pro¬ 
duction occurs predominantly around points where the axion velocity vanishes. For chaotic 
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Figure 6. The effective wavenumber k rescaled by the coupling C/f. The two black lines correspond 
to the maximum excited wavenumbers for A = — 1 (lower curve) and A = +1 (upper curve). The blue 
line corresponds to k = 0 and the green-dashed line corresponds to an intermediate curve for A = — 1. 
The red-dashed and brown-dot-dashed lines correspond to modes for A = 1 that get excited for the first 
time during and after inflation respectively. Particle production occurs whenever k crosses zero. Note 
that these particle production events occur in pairs which are closely spaced for high wavenumbers. 


inflation with a quadratic potential, </>(t*) ~ 1/3 during the first instance of </>(t*) = 0 and a 
simple application of Eqn. (5.13) and Eqn. (5.7) leads to 

n k « exp (-> fe < fc max~0.7y. (5.17) 

In the region of C/f 1 the particle spectrum approaches a step function in wavenumber- 
space with a sharp cut-off at k = hence the above formula, which suppresses the 

k— dependence, provides a useful estimate. Despite being a crude approximation, it provides 
an increasingly good estimate of the particle number for decreasing values of the combination 
m 2 /(C/f ) <C 1 and it can provide an upper limit on the particle production efficiency. 

We now consider the positive helicity state which becomes excited both during and after 
inflation. The range of modes excited during inflation is k < 0.7 C/f for chaotic inflation, as 
shown in the previous section, while the total range of modes excited is k < fc^ ax ~ 0.9 C/f, 
obtained from Eqn. (5.16). The difference of fc+ ax and ax plays a key role in our discussion 
of asymmetry in Section 5.2.3. 

An important characteristic of fermion preheating is the evolution of the Fermi sphere, 
defined in this case as the range of comoving wavenumbers that are (partially) filled with 
fermions as the universe expands during preheating. We again consider only the large coupling 
regime C/f^$> 1, in order for the WKB approximation to be valid. For a matter dominated 
universe, arising from an axion potential that is predominantly quadratic during preheating, 
the inflaton amplitude decays as (f> ~ a -3 / 2 and we may approximate 


4>o 


sin(f + 50) 
a 3 / 2 ' 


(5.18) 
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where 59 is an arbitrary constant phase. We then find a simple expression for the maximum 
excited wavenumber during a given production event 

1 C 

^max= W=T^O- (5.19) 

V « / 

Since the scale factor increases, the highest-momentum state that can be excited decreases as 
time progresses. For a locally quadratic minimum, the derivative of the effective wavenumber 
is given by 


dk 

1 

kX C , 

dt 
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and thus the particle number can be written 


n k ~ exp 



m 2 \ 

" WfWJ ' 


(5.20) 


(5.21) 


Note that the particle production rate decreases as a function of time. When taken together 
with Eqn. (5.19), Eqn. (5.21) implies the first few production events are the most important 
for determining the characteristics of the particle spectrum produced in this model. These 
events produce the most particles out to the largest wavenumber with the greatest efficiency. 
Subsequent production events can only affect smaller and smaller wavenumbers and with a 
decreasing strength. These subsequent events simply introduce a band structure that can be 
averaged out for all practical purposes, analogously to the static universe case. 

Finally, we note the difference between the case considered here and the case of Yukawa- 
coupled light fermions (m <C m§). In the limit of large Yukawa coupling, the maximum 
comoving wavenumber that is excited grows with the expansion as k mSLX ~ a 1 / 4 where a{t ) is 
the scale-factor, normalized as a = 1 at the end of inflation [14]. This result can be derived 
by noting that the particle production probability is given by 


n k ~ exp 



(■fc/a) 2 A 

hcfro/a 3 / 2 ) ’ 


(5.22) 


and thus the time dependence of the maximum wavenumber excited can be read off as k max ~ 
a 1 / 4 . In the case of heavy fermions (m > m^) the result is more complicated, giving the final 
form of n k as a shifted Gaussian in k (when studied numerically, this Gaussian is found to 
contain even more structure), rather than the smooth ferrni sphere of the light fermion case. 

As time progresses, the successive production events at each point where k = 0 change 
the particle number with decreasing effect. This is illustrated in the upper panels of Fig. 8 
where we show the average particle occupation number of Eqn. (5.11) as a function of time. 
We have so far only considered particle production from the first time where k = 0. In 
order to accurately estimate the final number densities for each helicity state, we extend our 
analysis to successive events. 

The problem of several particle production events can be mapped to the well-known 
quantum-mechanical problem of a particle undergoing multiple scatterings. In the static 
universe the constructive or destructive interference creates a well-defined band structure. 
However, in the expanding universe case, successive phases can be thought of as random. 
The calculation is analogous to that in [15], and we simply quote the result here, referring 
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Figure 7. Particle number for m = 1, C/f = 20 and A = —1,1 (left and right respectively) in 
an expanding universe. The different curves show the particle number at t = 0, at the time when 
all fc-modes of the A = 1 state have been produced for the first time and after 1,2,3 zero-crossings 
respectively from bottom to top. Curves are shifted vertically for clarity. The vertical lines correspond 
to k = 0.7 C/f and k = 0.9 C/f. 


the reader to the original works for more details. After N successive production events the 
smoothed particle occupation number is given by 

i= 1 

This result shows that Eqn. (5.7) can be used as to characterize particle production at all 
times, despite being strictly true only for the first particle production event. We test its 
accuracy in the next section. 

Fig. 8 shows the average particle number as a function of the successive k zero-crossings 
for a fixed value of the coupling C/f = 50 and varying fermion mass. The first result is the 
fact that particle production is suppressed for larger values of the mass, and effectively shuts 
off for m 2 /(C/f) > 0.5, and m 2 /(C/f) > 0.3, equivalently m > 5 and m > 4, for the positive 
and negative helicity mode respectively. A similar behavior was seen in Fig. 5 for the static 
universe analysis. 

The second important feature of Fig. 8 is the decline of the final particle number for small 
values of the mass. This presents a major departure from the results of the static universe 
analysis. In the static universe approximation, the particle number oscillates between two 
extrema for small masses, as shown in Fig. 5. From Fig. 6, note that particle production 
events always come in pairs of similar strength, especially for wavenumbers near /c max . This 
is shown in Fig. 6 by the two pairs of brown and green dots for the A = +1 and A = 
— 1 states respectively, which represent successive production events. For lighter fermions, 
particle production is very efficient, and the particles that are created by the first event 
tend to be destroyed by the second event. In the static universe the production efficiency 
is constant, but in the expanding universe subsequent events are much less efficient, as seen 
from Eqn. (5.21). Furthermore, as show in Section 5.2.2, particles are produced over a smaller 
range of wavenumbers as the universe expands. Taken together, somewhat surprisingly, the 
net result is a low overall particle number for light fermions. We also demonstrate this result 
analytically. Applying Eqn. (5.23) and taking the first two particle production events to be 


(5.23) 
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very efficient, that is e wpi = 1 — ~ 1, we get 

nf = ei + e 2 - 2eie 2 « ei + e 2 <C 1 . (5.24) 

meaning that (on average) modes that start almost fully occupied become almost empty after 
the second production event. 

An important qualitative feature of the expanding universe case is the similarity between 

( 2 ) 

the particle spectra after two production events, n, . and the result after a very long time, 
n[°°. Successive production events become increasingly less efficient and affect a decreasing 
range of wavenumbers, therefore these first two production events give the dominant contri¬ 
bution to the final particle spectra. We have already seen in the case of light fermions that 
these two events lead to the creation of a large number of particles followed by the immedi¬ 
ate destruction of many of them. The upper panels of Fig. 8 show that particle number is 
effectively constant after the first two oscillations. This is important because it means that 
any helicity asymmetry will be generated and finalized almost immediately after the end of 
inflation. 

Production events can be always described in pairs, therefore simple intuition tells us 
that, on average, modes with n& < 0.5 initially tend to grow their population towards unity 
after the second k = 0 point, while on average, modes with > 0.5 initially tend to reduce 
their population (away from unity). At fixed coupling C/f, production of lower mass fermions 
is more efficient, and they are mostly destroyed by the double production event. In the large 
mass region on the other hand, the efficiency of each event is lower, but successive scatterings 
are more likely to be constructive. Taken together, this suggests that the first two production 
events will produce a maximum of (n*,) ~ 0.5. This is exactly what the blue and red lines in 
the lower panels of Fig. 8 show. 


5.2.3 Asymmetric helicity production 

From our discussion of the range of excited wavenumbers and the evolution of the Fermi 
sphere, we already see a first hint of helicity asymmetry, simply from the fact that the positive 
helicity mode is produced up to a larger wavenumber ax > A;“ ax . With the average particle 
number in Eqn. (5.11) and the total particle number density 


n = 


/ 


d 3 kn k , 


(5.25) 


the particle asymmetry between the two helicity states becomes 

r^max /*^rnax 

An = n + — n~ « 47r / k 2 (n£ — n \“) dk + 47r / k 2 n//dk 
Jo J & max 

= Y(fcmax) 3 < n + - n ~)[+ Y {(k+ ax ) 3 - (fe- ax ) 3 ) (n + )ll2 • (5-26) 

For rn^cj) 2 inflation, the first term scales as (A:”^) 3 ~ 0.7 3 (C//) 3 ~ 0.34 (C//) 3 , while the 

second term scales as ((A^ ax ) 3 — (fc ~ ax ) 3 ) ~ (0.9 3 — 0.7 3 ) (C/f) 3 ~ 0.39 (C/f) 3 . The two 
terms can therefore be comparable, depending on the structure and values of n^. Having 
given estimates of the range of excited wavenumbers, we can now estimate the particle number 
and study the cases where it is possible to generate large asymmetries in the different fermion 
helicities. 
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Figure 8. Upper panels: Average particle number after successive production events for the positive 
(left) and negative (right) helicity states for C/f = 50 and TOf er mion = 0.5,1,2, 3,4 (blue, green, red, 
black and brown respectively). The 0’th production event occurs during inflation. Lower panels: 
Average particle number after the first (blue) and second (red) production event and after multiple k 
zero-crossings (green). 


Returning to Fig. 8, we note that positive and negative helicity states have compara¬ 
ble average particle occupation numbers. However, when converting this to a total particle 
number, one has to take into account the fact that each helicity state occupies a different 
volume in phase space. For m^cj) 2 inflation, for example, &;+ ax = 0.9 C/f for the positive he- 
licity state and fc“ ax = 0.7 C/f for the negative helicity state, which means that the particle 
number of the positive helicity state is more than twice that of the negative helicity state 
for ( n + ) = (n~). Thus the positive helicity mode is dominant, and a helicity asymmetry is 
present. We will quantify this asymmetry further using numerical simulations. 

5.2.4 Numerical results 

To test our analytic and semi-analytic results, we numerically evolve the equations of motion 
for the fermions in the homogeneous axion and gravitational held. We evaluate the evolution 
of the fermion number for a range of couplings 2 < Mp\C // < 500, and fermion masses, 
such that 0 < m 2 /(C/f) < 1. We neglect the back-reaction of the produced fermions on the 
axion evolution; back-reaction can be safely ignored for values of the coupling that satisfy 
(C/f)Mpi < 10 3 , at least during inflation and for the first stages of preheating, where ah 
interesting effects are expected to occur. 

Later axion oscillations (as shown in Fig. 8) become increasingly irrelevant, and so we 
focus on the first couple of fermion production events, which determine the final result to a 
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good approximation. For the first two production events Eqn. (5.23) becomes for the two 
helicity states 

2e ~ Mf ) (l - . 

"* ” ^ \ (‘ “ 2e “' W>a ) l 1 * . ( 5 - 27 ) 

which provide an increasingly good approximation for the numerically calculated result for 
C/f » 1, both after the second production event as well as after multiple ones. 

The important feature of the WKB-based results is their dependence solely on the com¬ 
bination m 2 /(C //), rather than on the fermion mass and coupling separately. However, as 
shown in Fig. 9, this is not strictly true for the particle number obtained from a numerical 
evolution of the system. There is a significant departure from the WKB result that occurs for 
wavenumbers that are near the maximum excited wavenumber. In particular, smaller cou¬ 
plings tend to have a larger “tail” in the nk distribution. This is the same behavior manifest 
in the static universe case. Because there is much more phase space at high wave number, 
the volume factor d 5 k weighs these modes more, and this discrepancy widens when the to¬ 
tal particle number density or average particle occupation number (n^) is evaluated. This 
phenomenon occurs for the positive and negative helicity states alike, leading to an enhanced 
(nfc) at smaller values of the couplings C/f, for the same ratio m 2 /(C/f). 

For order unity values of the coupling, the departure from the semi-analytic result is 
even larger - this is perhaps not surprising as the WKB approximation breaks down in this 
region. After several production events, these smaller couplings exhibit a slow but steady 
parametric excitation of larger wavenumbers, reaching k max ~ 2 C/f. Despite the occupation 
number being small njCl in this range, the particle number is again enhanced due to the 
d 5 k phase-space factor. Because we calculate the total particle density by integrating nk up to 
a sufficiently large wavenumber in order to enclose all particle production, and normalize it by 
^’max (given by the WKB approximation), the average particle number density {n±) can exceed 
unity for late times and small couplings. Furthermore, because of the steady excitation of 
these short-wavelength modes, the particle number keeps rising as time progresses, and there 
is no well-defined asymptotic particle number for late times, contrary to the case of larger 
couplings. A further effect of the behavior of these modes is the restoration of the symmetry 
between negative and positive helicity states. This leads to An = n + —n~ taking both positive 
and negative values for C/f = 2 in our simulations. However, the range of C/f = 0(1) is 
not of particular interest. This is due to the fact that the total number of particles (and the 
corresponding helicity asymmetry) is proportional to the cube of the maximum wave number 
^’max ^ (C/f) 3 , which implies that for a considerable effect to be achieved, we need C/f 3> 1. 
This makes model-building more robust, since A n/(C/ f) 5 is largely insensitive to the value 
of the coupling in this region. 

Fig. 10 shows the total particle number density asymmetry in helicity states An. This 
quantity scales with the cube of the coupling, so we have divided out by a factor of (C/f) 5 to 
put all curves on the same axes. The main results of this calculation are the dependence of 
the normalized helicity asymmetry An/ (C/f) 5 on the number of axion oscillations, as well as 
on the coupling C/f and the fermion mass m, or equivalently on one of these two parameters 
and the ratio m 2 /(C //). Fig. 10 illustrates that, for large values of the coupling C/f 1, 
the normalized helicity asymmetry is largely insensitive to the value of the coupling and is 
only determined by the ratio m 2 /(C //). Furthermore, there is no significant evolution of the 
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normalized helicity asymmetry after the second production event. We omit the case C/ f = 2 
from the right panel of Fig. 10, due to the absence of a well defined late-time asymmetry, as 
explained above. 

From a model-building perspective, for a fixed fermion mass m, the resulting asymmetry 
can be easily read from Fig. 10 for different values of the coupling. The maximum value of 
the asymmetry occurs for m 2 /(C //) ~ 0.07 and the full width at half maximum is given by 

o 

777 , 

0-02 < — <0.2, (5.28) 

which serves as a condition for the existence of a strong helicity-asymmetry. In order to 
provide a rough, but intuitive estimate of the maximum asymmetry strength, we start with 
the observation that the average occupation number for each helicity state, when considered 
as a function of m 2 /(C //), peaks at ( n fc) max ~ 0.5. The total particle number is 

n± = J i 3 kn t = („±) ^(*4(5.29) 


where we used the peak value (n) ~ 0.5. By assuming that the two distributions for the n/, 
and nV peak at the same value of m 2 /(C/f ) the asymmetry is simply 
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(5.30) 


Thus, for the peak value of the asymmetry, we see that the first term in Eqn. (5.26) vanishes, 
while the particle number density average in the second term can be taken as approximately 
0.5. For chaotic inflation with a quadratic potential, the values of lead to 


An m 2^,2 



(5.31) 


This formula provides a slight underestimation of the peak helicity asymmetry by about 10%, 
as seen in Fig. 10. However, both the accuracy, as well as the simplicity of this formula make 
it very useful for exploring the fermion production capability of different inflationary models, 
like the axion monodromy potential which will be discussed in the next section. 


6 Axion Monodromy Inflation 

In this section, we extend the results derived above in Section 5 for quadratic chaotic inflation 
to the simplest model of axion monodromy inflation with a potential given by Eqn. (2.3). All 
the machinery we have developed in Section 5 can be transferred directly to this case. The 
only differences are in the specific numerical values of \dk/dt\ at the various production events 
and the range of excited wavenumbers. 

We begin by rescaling our parameters to rewrite the background inflaton equation of mo¬ 
tion in a dimensionless form. To generate a spectrum of curvature fluctuations with an ampli¬ 
tude that matches the observed microwave background, we choose log 10 (/r/Mpi) = —3.22 [40]. 
Measuring the inflaton held in units of the Planck mass, and time in units of 1/y/Mp\/Jfi, 
the background equation of motion for the inflaton becomes 

(/ + 3H<j) 4- ^ = 0. (6.1) 

+ Vc 
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Figure 9. Particle number density after two production events for the negative (left) and positive 
(right) helicity state for chaotic inflation. Different colors correspond to different values of the axion- 
fermion coupling C/ f — 2, 5,10, 20, 50,100, 500 from top to bottom. The black-dotted line corresponds 
to the WKB result given by Eqn. 5.27. 
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Figure 10. Measure of the total helicity asymmetry An = n < - + ' ) — ' rescaled by the axion- 
fermion coupling cubed (C//) 3 after the second axion crossing (left) and after many (10) axion 
crossings (right) for chaotic inflation. Different colors correspond to different values of the axion- 
fermion coupling, consistent with Fig. 9, while the black-dotted line corresponds to the WKB result 
of Fig. 9. 


In what follows, wavenumbers will be measured in terms of yj 

Most of the qualitative features of the rn^cj) 2 case can be directly applied to monodromy 
inflation. For monodomy inflation, gravitational particle production during inflation is biased 
by the rolling axion so that only a single helicity is produced. After inflation both helicity 
states are excited at the points where k(t ) = 0, which differ for each helicity state and 
wavenumber. The range of excited wavenumbers is smaller for each subsequent production 
event, due to the decay of the inflaton oscillations, making the first two production events the 
dominant ones. For axion monodromy inflation, there is an additional parameter, <f> c , which 
controls the shape of the potential near the end of inflation and quantitatively changes the 
results. 

For the first few oscillations dominate and define the range of the Fermi sphere, 

as we found in the rn^cj) 2 case. The corresponding maximum wavenumbers excited during the 
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first production event for each of the two helicity states are shown in Fig. 11, along with 

/»;, \ 3 (^max) max) (n 

^AA/vmaxJ — (C //) 3 ’ v D-Z / 

which tracks the difference in phase space volume that is populated in each case. Note that for 
m 2 (j) 2 inflation (Afc max ) 3 ~ 0.4, which is similar to the values for axion monodromy inflation, 
in the range of 0.1 < (j) c < 1. 

For large values of cj) c >, the minimum of the potential behaves more and more like 
a quadratic model during the oscillation phase, and we find similar quantitative results to 
inflation. Specifically for (f> c = 1 the form of a(t)<p(t) is very similar for quadratic and 
axion monodromy inflation, with a difference in amplitude of about 25%. This leads to a 
difference in the number density of produced fermions of about 1.25 3 ~ 2. The maximum 
wavenumber excited in each axion oscillation, /q^ ax , grows with decreasing (j > c , as shown in 
Fig. 11. However, the combination (/c+ ax ) :! — (fc“ ax ) 3 decreases for <j) c < 0.3. For small values 
of (f c , both helicities are produced in increasingly equal amounts, resulting in smaller helicity 
asymmetries. The reason for this behavior can be seen from the (unphysical) limit (j) c —> 0. 
In this limit, </ >(t)a(t ) oscillates with a constant amplitude generating a Fermi sphere with 
a constant radius in momentum space. Each axion oscillation therefore produces particles 
out to the same maximum wavenumber with an efficiency that remains essentially constant. 
Hence for cj) c < 0.1 we have an increased particle production (due to the increased range of 
produced wavenumbers) and a decreased level of helicity asymmetry. 



Figure 11. Left panel: The maximum initially excited wavenumber for A = +1 (upper blue curve) 
and A = — 1 (lower red curve) is shown in units of C/f for axion-monodromy inflation. Right panel: 

Plotted is the ratio (fc+ ax ) — (l max ) 3 / (C/f) 3 , which can be used to estimate the helicity asymmetry 
from Eqn. (5.30). 


The average particle number after the second production event for each helicity state is 
shown in Fig. 12. This is a very good approximation for the final particle number, especially 
near the peak of the curve (nT) as a function of m? j(C //), where both the maximum particle 
production and the maximum asymmetry occur. We present results for the range 0.2 < <f> c < 
0.6, where the asymmetry is largest, as expected from Fig. 11, but simulations in the whole 
range 0 C E [0.1,1] show similar results. As expected, the maximum average occupation 
number is (nA) ~ 0.5. Making use of Eqn. (5.30) and Fig. 11, for axion monodromy inflation, 
the total helicity asymmetry in the number density of left- versus right-helicity fermions is 
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given approximately by 


An 


max 



(6.3) 


where the subscript refers to the maximum asymmetry. The value of m 2 /(C/f ) for which 
the maximum asymmetry occurs is a model-dependent quantity, but for the model of axion 
monodromy inflation studied here m 2 /(C/ f) ~ 0.05. The result of Eqn. (6.3) is valid, with 
good accuracy, for 0.1 < <p c < 1. The helicity asymmetry in the monodromy case behaves in 
a very similar fashion to that in the m 2 cj) 2 case, however, with a different magnitude and a 
different value of the ratio m 2 /(C/f) for which the maximum asymmetry is achieved. Both 
the exact magnitude of the maximum asymmetry and the corresponding value of m 2 /(C/f) 
depend on (f> c . 
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Figure 12. Average particle number after two production events for the case of axion monodromy 
inflation for <f> c = 0.2, 0.3, 0.4,0.5,0.6 (blue, red, green, brown and black respectively). The maximum 
value of ( n ± ) « 0.5 is evident. 


7 Conclusion 

In this work we have studied the behavior of derivatively-coupled, massive fermionic degrees 
of freedom during, and immediately following axion-driven inflation. When the masses of the 
fermions are degenerate, these fermions pair up to make Dirac fermions and the coupling is 
to the axial-vector current. 

During inflation, the motion of the inflaton is monotonic and the coupling of the fermions 
to the rolling axion acts as an effective chemical-potential for helicity. This chemical potential 
biases the gravitational production of helicity states of the fermions. In the absence of the 
coupling, massive fermions are produced gravitationally by the expansion of spacetime equally 
in all helicity states. This production is most efficient for light fermions; the production rate is 
exponentially suppressed by the ratio of the fermion mass to the Hubble rate. When the axion 
coupling is turned on, we have shown that one helicity state is produced more efficiently while 
the other is strongly suppressed. This mechanism allows for the gravitational production of 
heavy fermion states that would otherwise be heavily suppressed. 
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Following inflation, as the axion oscillates, the effective frequency of the fermion field 
varies non-adiabatically resulting in particle production. These oscillations also mean the ax¬ 
ion velocity changes sign, resulting in the production of both fermion helicities. The damping 
of the inflaton amplitude after the end of inflation, due to Hubble friction, results in biased 
fermion production in favor of the helicity state that is excited during inflation. This is due to 
the declining efficiency of the production rate as the axion damps combined with the fact that 
production of the second helicity does not begin until the axion velocity changes sign. For 
certain combination of the fermion mass and axion-fermion coupling, the produced fermions 
can have a significant helicity asymmetry. This helicity asymmetry is robust to the details of 
the reheating history, at least in the limited models we considered here, because it requires 
the amplitude of the axion oscillations to decay faster than l/a(t), which does not need any 
significant tuning of the model. 

We studied the post-inflationary evolution of the particle number in both quadratic 
chaotic inflation - rn^cj) 2 inflation - and in the simplest model of axion monodromy infla¬ 
tion. While there are minor numerical differences between the models, broadly we find very 
similar behavior. In both cases the average particle number is solely defined by the com¬ 
bination m 2 /{C //) in the regime of large coupling, while there is increasing deviation from 
this behavior as we decrease the coupling. Further, in both models the range of excited 
wavenumber shrinks with time as well as the production efficiency of any single event, leading 
to the first production events giving the dominant contribution to the final particle spectrum. 
The average particle number remains almost constant after the second production event for 
C/ f 1. Finally, the maximum achievable helicity asymmetry in the number density of pro¬ 
duced fermions, An, scales as (C//) 3 , with a proportionality factor of 0(1), which depends 
on the specific inflationary model. 

Neutrino helicity is equal to lepton number in the standard model of particle physics, 
therefore the asymmetric production of helicity may be important for the generation of the 
matter-antimatter asymmetry in the Universe. We apply the results of this work in a compan¬ 
ion paper [12], where we propose that an axion-inflaton coupled to left-handed standard-model 
neutrinos as a mechanism for the observed baryon asymmetry in the Universe via leptogenesis. 
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A Some 2-component spinor notations 

In this Appendix, we briefly review the 2-component spinor notations and conventions we have 
used in this paper. Our notations follow [35], to which we refer the reader for more detail. We 
will use un-dotted indices to denote fields transforming under the left-handed representation 
of the Lorentz group, (^,0) (ip a —> and dotted indices to denote fields that transform 

under the right-handed representation of the Lorentz group, (0, |), ('i/’t —> ( M )• There 
are two additional irreducible representations of the Lorentz group which are equivalent to 
the left- and right-handed representations above. The fields that transform under these 
representations are denoted by raised indices, i/j 01 and ifj a . It is convenient to consider ?/>" as 
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a row vector, and ip a as a column vector, while ift 01 = (V ,a )^ as a column vector and tjy- is 
a row vector. We will use the summation convention that repeated pairs of upper and lower 
indices are summed. 

These spinor indices are raised and lowered using the spinor metric tensors, denoted by 
the 2D antisymmetric e symbol, where our conventions will be 

e 12 = —e 21 = €21 = — ei 2 = 1, (A.l) 


and we formally define = (e aj g)* and e a @ = (e Q ^)*. Note that these conventions imply 
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€ — 6 € ^a/3i ^7 S — ^7 a^S/3^ 

These epsilon symbols satisfy 

^=-pa+s‘ji =- s H + s H’ 

and 

f nf ^ 7 ~ P € ' f ^ 7 — P 

€a/3€ — 0 a , — 0 a . 

We use the sigma, or Pauli matrices, cdT and a aot °, defined by 


(A.2) 


(A.3) 


(A.4) 


< 7 ° = < 7 ° = 


1 0 
0 1 


2 -2 
cr = —cr = 


0 -* 

! 0 I ’ 


1 -1 
cr = —cr = 


_3 -3 

cr = — cr = 


0 1 
1 0 

1 0 
0 -1 


(A.5) 

(A.6) 


These sigma matrices are hermitian, and are defined with a contra-variant (upper) spacetime 
index. Denoting the 3-vector of Pauli matrices a = (a 1 , a 2 , a 3 ), the Eqn. (A.5) is equivalent 
to 


cd‘ = ( 12 x 2 ,*?), ^ = (12x2,-07 (A. 7 ) 

For convenience, we work only with left-handed fields, that is, fields that transform 
under the (|,0) representation of the Lorentz group. Since these representations are related 
by hermitian conjugation, right-handed spinors are then simply conjugates of left-handed 
spinors. 
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